library(dplyr)
library(ggplot2)
library(stargazer)
library(geodata)
library(terra)  
library(sf)
library(stringr)


load("data_analysis.RData")

## Keep only if case entered during Santos
base_foranalysis <- base_foranalysis_april24 %>%
  filter(min_fecha_adminis < as.Date("2018-08-07"))

base_foranalysis$days_until_exit <- ifelse(
  is.na(base_foranalysis$tiempo_admin_salieron_registro),
  base_foranalysis$tiempo_admin_if_trapped,
  base_foranalysis$tiempo_admin_salieron_radi
)

base_foranalysis$trapped <- ifelse(
  is.na(base_foranalysis$tiempo_admin_salieron_registro),
  1,
  0
)

base_foranalysis$tiempo_admin_both <- ifelse(
  !is.na(base_foranalysis$tiempo_admin_salieron_radi),
  base_foranalysis$tiempo_admin_salieron_radi,
  base_foranalysis$tiempo_admin_if_trapped
)

base_foranalysis <- base_foranalysis %>%
  mutate(depto_name = toupper(trimws(departamento)))   

colombia_dept <- gadm(country = "COL", level = 1, path = ".")

colombia_dept_sf <- st_as_sf(colombia_dept)

colombia_dept_sf <- colombia_dept_sf %>%
  mutate(depto_name = toupper(trimws(NAME_1)))


dept_stats <- base_foranalysis %>%
  group_by(depto_name) %>%
  summarise(
    n_cases       = n(),                                 
    mining_titles = sum(num_titulos, na.rm = TRUE),      
    .groups = "drop"
  )


dept_data <- colombia_dept_sf %>%
  left_join(dept_stats, by = "depto_name") %>%
  mutate(
    n_cases       = ifelse(is.na(n_cases), 0, n_cases),
    mining_titles = ifelse(is.na(mining_titles), 0, mining_titles)
  ) %>%
  filter(!grepl("SAN ANDR", depto_name))   


dept_centroids <- dept_data %>%
  st_point_on_surface() %>%
  dplyr::select(depto_name, n_cases, mining_titles, geometry)

offices_depto <- toupper(c(
  "Antioquia",
  "Bolívar",
  "Caquetá",
  "Cauca",
  "Cesar",
  "Chocó",
  "Córdoba",
  "Magdalena",
  "Meta",
  "Nariño",
  "Norte de Santander",
  "Putumayo",
  "Tolima",
  "Valle del Cauca",
  "Bogotá D.C."
))

label_points <- dept_centroids %>%
  filter(depto_name %in% offices_depto) %>%
  mutate(
    lon = st_coordinates(.)[, 1],
    lat = st_coordinates(.)[, 2]
  ) %>%
  mutate(
    # X offsets
    lon_label = case_when(
      depto_name == "ANTIOQUIA"           ~ lon - 0.8,
      depto_name == "VALLE DEL CAUCA"     ~ lon - 0.8,
      depto_name == "CAUCA"               ~ lon - 0.6,
      depto_name == "NARIÑO"              ~ lon - 0.6,
      depto_name == "CHOCÓ"               ~ lon - 0.6,
      depto_name == "CÓRDOBA"             ~ lon - 0.6,
      depto_name == "BOLIVAR"             ~ lon + 0.8,
      depto_name == "MAGDALENA"           ~ lon + 0.6,
      depto_name == "CESAR"               ~ lon + 0.8,
      depto_name == "NORTE DE SANTANDER"  ~ lon + 0.8,
      depto_name == "META"                ~ lon + 0.6,
      depto_name == "PUTUMAYO"            ~ lon + 0.6,
      depto_name == "TOLIMA"              ~ lon + 0.4,
      depto_name == "CAQUETA"             ~ lon + 0.6,
      depto_name == "BOGOTÁ D.C."         ~ lon + 0.4,  
      TRUE                                ~ lon
    ),
    lat_label = case_when(
      depto_name == "CHOCÓ"               ~ lat + 0.5,
      depto_name == "MAGDALENA"           ~ lat + 0.4,
      depto_name == "BOLIVAR"             ~ lat + 0.4,
      depto_name == "NORTE DE SANTANDER"  ~ lat + 0.4,
      depto_name == "META"                ~ lat - 0.3,
      depto_name == "CAQUETA"             ~ lat - 0.4,
      depto_name == "PUTUMAYO"            ~ lat - 0.4,
      depto_name == "NARIÑO"              ~ lat - 0.4,
      depto_name == "TOLIMA"              ~ lat + 0.2,
      TRUE                                ~ lat
    )
  )

p_cases <- ggplot() +
  geom_sf(data = dept_data,
          aes(fill = n_cases),
          color = "white", size = 0.2) +
  geom_text(
    data = label_points,
    aes(x = lon_label, y = lat_label, label = depto_name),
    size = 3,
    check_overlap = TRUE
  ) +
  scale_fill_gradient(
    low  = "grey90",
    high = "darkred",
    name = "Restitution\ncases"
  ) +
  coord_sf() +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    axis.text        = element_blank(),
    axis.ticks       = element_blank(),
    legend.position  = "bottom",
    legend.box       = "horizontal"
  ) +
  labs(
    x     = NULL,
    y     = NULL,
    title = "Restitution cases by department"
  )

p_cases


p_mining <- ggplot() +
  geom_sf(data = dept_data,
          aes(fill = log(mining_titles + 1)),
          color = "white", size = 0.2) +
  scale_fill_gradient(
    low  = "grey90",
    high = "darkblue",
    name = "Mining titles (log)"
  ) +
  coord_sf() +
  theme_minimal() +
  theme(
    panel.grid.major = element_blank(),
    axis.text        = element_blank(),
    axis.ticks       = element_blank(),
    legend.position  = "bottom",
    legend.box       = "horizontal"
  ) +
  labs(
    x     = NULL,
    y     = NULL,
    title = "Mining titles by department"
  )

p_mining



ggsave(
  filename = "map_restitution_cases.pdf",
  plot     = p_cases,
  width    = 6,
  height   = 8,
  units    = "in"
)

ggsave(
  filename = "map_mining_titles.pdf",
  plot     = p_mining,
  width    = 6,
  height   = 8,
  units    = "in"
)
